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Abstract 


This paper explores an idea of S Kant nf 

attitude of a satellite, given that the cravitv fi m * d geometer measure the 

measure a eo„bi«to" f ^ “d^ -T ‘"T 7 S ”<* ^adiometem actmdly 

from obvious. The paper denuSSS ™ ZTfj *? •—». the anmter ia f„ 

is dynamic estimation, based on the momentum biased°P t accl ^ racy ' The technique employed 
nominally pi*., p„ inlrf| „ K , , , J™'"'™ "““O' 1 Eul« equatmna. The satellite is assumed 

The attitude estimator is unusual While the standi d ®" d partly r& ndom drag torques, 

is t~d, the feedback gain m^mU. 

a measure of the terminal covarimii of the^rmTn 5 7 lh “ Iy ' If. chcen to minimise 
«d the po«er spectra of the p T f d ' P '" da «“ 
multiple solutions of Lyapunov equations. ‘ 1Ses - An lnt *gration is required over 

1 Notation &; Units 

Uppercase bold roman letters are 2 dimensional arrays- e e F r„„„ u m 

are column vectotu; e g ., r . Magnitnde3 „ f vect0[s bold r»m.n or greek letters 

are indices. The Einstein summation convention is used f ’ , g j F ' lowercase greek subscripts 
signify time derivatives; e.g., £ = dx/dt \ T for repeated lower case greek indices. Overdots 

variables. Sines and cosines are denoted by s and c respeTuvely^ 0168 tranSP ° Se ' Primes denote sc ^ 
A = constant matrix in Riccati equation 

fy^-T eXte , rnal non ~gravitational acceleration on spacecraft 
a, - inertial acceleration of the ith accelerometer ? 
a - process noise state distribution matrix 
^ vector of state concern values* C — +’ 

DM = matrix satisfying L y, pun<> ; 

F = plant rrmtrbcf ^ K py ?' J “ , UIUt VeCt0r along axis Q in coordinate system x 
C = „1T 7 ’ fe - external non-gravitational force on spacecraft 

_ = OT ® rea gravitational constant = 6.67259 x KT 11 N-m 2 /ke 2 

K = filter feedh p • der ~ 0Ve £, aI1 spacecraft inertia tensor 

= coiners kr- H Jr M d' 1H f in Rieca,i e ‘ ,u * ,i ° n 

M = combined -equivalent- ‘J'te £.t2X d^ta ^ ™ 

m = spacecraft mass; also field source mass ' 

Q(«) » funcrionofpler defi f m ° f the 681101816 

firfieM" 111000 ^ 1 ^ 10 " matriX WUh deky ^ W) = averl g °e r power ** Pr6SSUre 

S(ut) = gen 0 e S rai 0 noLe poVefsp^trum- S ^ ° f pressure 

s, t superscripts signify spacecraft ank tmiecloiyly^rj^ 

J = time m seconds; f' = t/C t = scaled time; t 3 Z filter sett lT ne tim. 

V = BS C U^- 0IS w e h r aSUrement diStribution matrix : « = vector of controls 
~ whlte Process noise effect matrix 

oo ' «-“»>«■ 
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v 0 = satellite orbital speed 

w x* _ KU = process noise effect matrix 

pmc^LJvector; = dimensionless air drag random process 
x = F - VM-'H = linear term matrix in Riccati equation 
x = state vector; x = estimate of x, x — dx/dt 

Y = measurement noise distribution matrix measurera ents 

2 = F - KH = observer system matrix, z - veci 

ttZZSSZX r 0 - 0 re/c 3 - gradient scalar due to mass - a. distance r 

A(Z) = eigenvalue of Z; » - *W P "‘ ° h “SsM X 10“ m 3 /s 3 

C = Om. = gravitational non _p, v itati„nJ external torque 

< . X - x = error tn the « » f[equ S,cy used in power spectra 

(jj = spacecraft angular velocity, g 11 , , , me an motion 

w = break frequency in power spectrum; u> 0 = orbital 

qt xj nwPV pr we have also followed common 
Unless otherwise stated, the units used .in , t is gradient. The natural SI unit is (m/s 3 )/m, or 

practice in the held of gr.diome.ry on Ih ‘ “ „„ lhe ordet of 1.5 X 10- s‘ 3 , and are 

iust s -2 Since gradient components at the earth ... There has now been world wide 

^tU measured to 10- s- or b-er this has proved " e Everywhere. n the formulas; 
acceptance of the Ebtvos unit: 1 E = 10 s • Here 
but Ebtvos units will be generally employed in the text. 

2 Static Attitude Estimation 

The gravitational potential due to a particle of mass m at a distance r is: 

$ = -Gmjr ' 

The vector gravitational held at this point, due re re, is the acceleration of a free test particle there: 

g = _V4> = —Gmr 3 r s -lor 

Finally, the gravity gradient tensor field due to m, is. 

r = v g = r 0 (yr - l3 ) 


( 3 ) 


w 2 =ro = Pe /r J 


( 4 ) 


a r .v. nnriVi The actual potential of the earth is complicated; 
in which p e is the gravitational ‘ :ons ^" 0 ® - n low earl h orbit, less at higher altitudes. The variations 

but differs from (1) by only about 1 par Thus if spaC ecraft attitude is actually inferred from 

in turn are known to better than 1 par . of ’ the field wou ld lead to corresponding attitude 

gradiometer measurements, this error mk ^ not the worst error contribution. In any 

^Tl the accuracy with which a gradiometer can measure attitude, ,.«en 

C faat the field is known] so the study neglects hold knowledge errors. 

On the other hand, neglect of the known intent is to determine 

the spaceerati ° rbit is isk "' here as circu,M . h 
^ Kpd hv sets of right handed orthonormal base vectors e Q , where 

«• - — ti systOT Thl> “ the 
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SK? Z71VctTZt 10 r' mmeKr inpu * and •» »“>» >«"«=. ™ 

gned. for simplicity, it will be assumed that the origin of e a is at the spacecraft center of m «=c Th« 

e 2 completes a right handed system, and is along the spacecraft velocity vector e‘ rotates uniformlv »t 
. me aboui .$ relative to inertial system that won't need to be fden.ified funb.r J 


vw vv lutiliuicu 1 LU tit 

The connection between systems may be described by a matrix of direction cosines A: 

e a = Aapep 

r-= 


(5) 


A = 


1 9 -<fi 

—9 1 ip 

<t> -i> l 


The need for e is that the earth fields g and T are most conveniently expressed there: 

g‘ = -r 0 re‘ = r 0 r[— 1, 0, 0] T 
r f =r 0 diag[ 2 , - 1 , - 1 ] 

and expressing these in e*, where the instruments reside: 


s’ = Ag* = r 0 r[-l, 9, 


r°= Ar t A T =r 0 


2 

-39 

3<p 


-39 

-1 

0 


3 <p 

0 

-1 


( 6 ) 


(7) 

( 8 ) 

(9) 

( 10 ) 


dlTnt 0 , PiTtVv thl t,”h leS N ° le fct ' ,h “' 7 hile PitCh “ d ™" turn “!> fa expressions, yaw 
aoes not. Physically, this is because r is an axis of symmetry of the fields. 

If we could measure either g or T in e*, we could infer both 9 and 4 >. Alas accelerometers don’t 
measurement of FJ 3 could somehow be made, the error in 4 > would be: P ’ P 


( 11 ) 


— SF/ (3Fo) + 3 4>6r/r 

Suppose an orbit altitude of 500 km. Then r = 6 867 x 10 6 m and r„ - ion v \ .. „ 

of E r. wo te ih h “" in ;s u, 5 e r 8 * i<rs rad “* ,ysis •' 

of „ r is- r " each th - 2 "'> contribution to the error comes from the uncertainty in the 
Since s fin Upposmg< ,r- 10 m, and <f = 0.1 rad, this contribution to £41 comes to 4 37 x 10" 7 rad 

•he traei^SJSr 

With an 0.5 m separat.on and independent errors, their required accuracy would be aCCelerometers ’ 

<5a = 0.5(10" n )/2 1/2 = 3.536 x 10~ 12 m/s 2 

within the capability of the best room temperature accelerometere today, operating in space. 

3 Dynamic Attitude Estimation 

If grad, ometers actually measured the gradient, then a model would be something like z = r plus noise 

in the e bS . °! r?J nP a° nent j- A leaSt SqUares ana 'y sis would ^en yield the covariance oHhe errors 
m the estimate of the d f ired * and 9, for each discrete sample z. However, as any red gradbmeL 

measurement z contains functions of u, and *, least squares analysis won’t suffice; and we have to resort 
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^ ^ shown ho. .he — «— 

of the errors in the estimates may be determined. A few results are given. 

!n°L ,f aligned .long the ej Supping . uniform den,.ty p, the sp.cecn.ft m« ». 

. . , ( 12 ) 

m = pi 1*2*3 

A typical density might be p = 1000 kg/m 3 ; and the principal moments of inertia are: 

Jj = m(l% + 13) / 12 ; J 2 = m(i? + /|)/12 5 J 3 = m(i? + i£)/ 12 

The orbit is assumed circular, at a radius r. Assuming an altitude of 500 1 km r = * 867 * , l f m ’ 
Wo = .0011095 rad/s, and T 0 = 1231 E. Also, the spacecraft speed in orbit v a o / • 

In [4] it’s shown that the Euler equations of rigid body motion, when modified to include an arbitrary 
bias momentum hw, can be written as: 


(13) 


Ju> = (Ju> + hw) Xu I + Tgg + Tf 


(14) 


^=^-S»te6r 

of the body derivatives, a much simpler procedure is to define the variation e by. 


OJ = U>o e 3 + 6 


(15) 


Another simplification comes by arguing that, in an earth pointing satellite, bias momentum, if any, is 

usually confined to the pitch axis: s (16) 

hyv = he 3 

Additional wheels for control aren’t precluded; it’s only required that their nominal momentum is zero. 
Substituting these relations into (14), and deleting quadratic terms m e, yields 

(17) 

(18) 

(19) 


Je = Wo (Je) x ej + u, 0 (Je£) x ( Wo e$ + c) + he’ 3 x e + r 9S + r e 

We also need r gg . The well known formula in e‘ may be put in the form: 

Tg g = 3r 0 ei X (JVj) 


Since only J> is readily available, and as what we really need is r-„, we need to work out 
T * = 3r 0 A[ei x (A T J*Ae t 1 )] = 3r 0 


(Jll — + ^23# - ^13 

[ {J\l - ^ 22 )^ + ^230 + J \2 J 


Note that, while nothing depends on J, there is a yaw t.r^ 

These also produce bias torques in roll and pitch. That s why, lor earth pointing 
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preferable to point some principal axis up. Moreover, by making this axis (ef) have the least J, r gg is 
restoring. Here, where the main issue is observability, it’s assumed that this condition is met, when 

J* = diagfj 1? J 2 , J 3 ] (20) 

In the examples, it’s further assumed that Ji < J 2 < J 3 , known to be the best configuration for gravity 
gradient stabilized satellites. With the principal axis assumption, the torque reduces to: 

r’ gg = 3r 0 [o, (j t - j 3 )t, (j, - j 2 )e f (2i) 

On putting this into (17), and expressing it in standard form, the Euler equations become: 

i\ — k\t 2 + 1 T e j 

h = k 2 (i + k 3 <f> + J 2 1 t k2 (22) 

€3 = k 4 6 + 1 r e3 

in which the constants are defined as: 

*1 = MJ 2 -^3 )-A]/Ji ; k 2 = MJ3- Ji )+ h]/ J 2 ; k 3 = 3r o (J!- J 3 )/ J 2 ; k 4 = 3r o (J x - J 2 )/ J 3 (23) 
If the gravity gradient and other external torques are neglected, then e, and t 2 decouple from £3 in (22), 

resulting in a harmonic oscillator with frequency w/v given by: 

w jv — ~k\k 2 (24) 

This is the natural nutation frequency, arising mainly from the momentum bias h. 

To complete the plant equations we must add the kinematical relations. With the same linearizing 
assumptions, these are easily shown to be (see for instance [4]): 

^ — £ i + u 0 4 > ; 4 , — e 2 - u 0 ip ; 0 = e 3 (25) 

We now have a linear system of plant equations of 6th order in c, %l>, <j>, and 6. 

The random process appearing in the Euler equations (22) is the external non-gravitational torque r e . 
At 500 km, this is largely due to air drag; and the random component is largely from variations in air 
density p a . For gradiometer studies, a flat earth barometric model was adopted in [4]: 

Pa(r + Sr) = p a (r)e~ 6r / h ‘ ( 26 ) 

where h 3 is the density scale height. At 500 km, [9] lists Pa = 1.905 x 10" 12 kg/m 3 , h„ = 83,000 m, and 

a mean free path of 25,000 m. These numbers are, admittedly, quite shaky. In any case, the dynamic 
pressure then comes from the speed: 

9 = p a vl/ 2 ( 27 ) 

and with the above numbers, 9 = 1.106 x 10~ 4 N/m 2 . Since the speed is along e', and the spacecraft 
attitude is not far from nominal, the steady force from air drag is very nearly: 

fe~ -qhl 3 C D e 2 (28) 

Because the mean free path is much larger than the spacecraft, drag is essentially Newtonian, with a 
coe cient Co — 2. However, since some inelastic, oblique, and diffuse scattering of air molecules is 
ley, t is Cd may be high, and Cp = 1.5 is adopted. We should also consider radiation pressure. 
Corresponding to q is I 3 /c, where I a = 1360 w/m 2 the mean insolance outside the earth, and c is the 
speed of light. Thus, the mean “radiation dynamic pressure” is 4.54 x 10 -6 N/m 2 , well below q\ and as 
the variations are much slower than for air drag, radiation pressure is ignored. [4] goes on to develop a 

statistics 1 model. It supposes that Pa is actually the mean of a distribution, to which a random component 
is added: 

Pr = PaWd{t) ( 29 ) 
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iua(t) is a dimensionless, zero mean, random function of position and time. At satellite speed, the spatial 
variation dominates. Suppose that w d (t) has a standard deviation <r w . Still, we need a power spectrum. 
Physically, we are looking at dynamic variations in density, with scale lengths of order h s , plus the 
orbital frequency variation due to solar heating of the atmosphere. The latter, while reaching substantial 
amplitudes, is confined to such low frequencies as to have little effect on the attitude estimates, and is 
ignored. As for dynamic variations, we can imagine variability on all length scales, but petering out 
below distances of order h t . This situation led to the development of the cubic power spectrum in [2]. 




(30) 


and zero otherwise. Suppose the autocorrelation of variations falls by half at a distance ah t . The time 
to travel this distance is A = ah,/v 0 , and [2] shows that, for the cubic spectrum, we should choose: 

(31) 


tdcw — ^ * 


2A 2ah„ 

We must also pick R w ( 0) and a. The best information presently available to us is an analysis of CACTUS 
data in [10]. Accelerometer data over approximately 800 s intervals was analyzed at altitudes between 
270 and 320 km. Density variations of ~ 4%, peak to peak were typical; rising sometimes to ~ 15%, 
during severe magnetic disturbances. The corresponding values are .014 and .05. A reasonable balance 
between these values would be ~ .02; but, allowing for a bit greater variability at higher altitudes, we 
have taken <r w = .025. Then, as these time series meet the oversampling conditions discussed in the 
Appendix R (0) = a 2 = 6.25 x 10~ 4 . As for a, [10] doesn’t show a power spectrum, but does give 
representative time series of a normal and a disturbed interval; and states that the apparent wavelengths 
concentrate in the range of 700 to 1500 km. Examination of the time series suggests that R(t) falls to 
0.5 at r ~ 50 s. Translating to our altitude, the corresponding distance is 381 km, when a = 4.6. Since 
for a sinusoid, R(t) falls by half at 1/3 of a wavelength, these numbers are at least consistent. Again, to 
allow for a bit more variability at 500 km, we have taken a = 4, leading to u cw = 0.03606 rad/s. 

It remains to convert this to torque. The overall drag force is very nearly. 

f e = -k f [l + w d (t)]el 


(32) 


kf = p a vlhhC D l2 

jposing a center of pressure at a location r cp in the spacecraft, the torque due this is: 

_ ..r_L.fi! ....MIL „ n — T- A T 


(33) 

(34) 


Note that there a deterministic bias force and torque, which must be treated correctly m the filter. Also, 
while our box structure has no torque along e * 2 , an actual spacecraft would likely have a small propeller 
torque on this axis. To allow for this below, a component r cp2 replaces the zero in (34). 

4 Measurement Model 

In [5] and [3], the instrument was modeled as measuring elements of the “intrinsic” tensor: 

T = r + w 2 I 3 - ww T + eu> (35) 


where e is the 3-index permutation symbol. The quadratic u> terms are centrifugal effects. Because 
the instrument is fixed in e 4 , there is no coriolis. Here, the instrument is dissolved into its component 
accelerometers, partly to avoid the noise correlations required in [5] and [3], but mainly to prepare for 
later studies. The gradiometer is taken as an array of 3 axis accelerometers, with input axes aligned 
along the e 4 . For entering symmetrical arrays, it’s convenient to identify a “center” of the instrument 
r c , relative to the origin of e 4 . Then, the tth accelerometer will have a position r oi , relative to the center. 
Thus, its location relative to the center of mass is: 
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For a perfectly circular orbit, the center of mass is subject to — uj^reJ. As for rotation effects, e* is 
rotating at a rate u>, relative to an inertial frame e". So, purely due to rotation, the inertial velocity of 
the ith accelerometer is (the superscripts indicate the frame in which the derivative is observed): 


d n d s 

T i — 3“r» = — r* + u x r, = w x r, 
dt eft 


the latter because r, is invariant in e 9 . Going to the next derivative 


. d n 

r * — — r< -WXTj+WX —Ti = U> X Ti + U> X (u> X Tj) 


(37) 


(38) 


Note that u> is the same, whether viewed from e" or e 9 . Finally, on including the external non- 
gravitational acceleration a c , the ith accelerometer is subject to: 

&i = -wlre\ + u> x (u> x r<) + cb x n + a* (39) 

On the other hand, the gravitational acceleration of the ith accelerometer is g* plus the correction at r, 
due to the gradient. From (7) and (4), this comes to: 

Si = -ulre\ + rTi (40) 

Actual accelerometers measure only non-gravitational acceleration; i.e., the difference between inertial 
and gravitational acceleration. These are identical in free fall, when an accelerometer measures zero. 
Conversely, an accelerometer on a table on earth measures the acceleration imposed by the table that 
keeps the instrument from falling through the floor. Thus, the ith accelerometer model is: 

z^ai-g^ + Vj (41) 

where v, is the noise in the 3 measurements. On substituting from above this becomes: 

z, = u> x (w x Ti) + tb X Tj — rTi + a* + Vj (42) 

Note that the acceleration of the center of mass has dropped out. The next step is to linearize this using 
(15). On neglecting the quadratic terms, and recalling that cb is the same in e n and e 9 , we get: 

Zj = w 0 (w c e3 + c) x (eg x r.) + w 0 e3 x (e x r^) + k x r* - Tr, + a,. + v* (43) 

We 11 work this out term by term, in the form of matrices of constants times the state variables, plus 
whatever is left over. Starting on the left: 

e 3 x ( e 3 x r ») = ^*363 — Ti = — [r**i , Ti 2, Of (44) 

(45) 


(46) 


e 3 x (e x r.) = r^e-egr^ 

The c term can’t be expressed directly in the state variables; however, from (22), there follows: 



1 

0 

0 

1 



= 

0 0 -r i2 




r a r i2 0 



. 

f r i3 0 -Til 1 




0 

w 

1 

jl 

to 


^2 

L 

0 0 0 J 
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^2 r i3 

0 

^3 r i3 

^4 r t2 




J 2 lf t3 r c2 

- J 3 1 r i2 r e 3 

e x Ti — 

0 


0 

k 4 rn 


^2 

<t> 

e 

+ 

*?3 r tl r c3 

- Ji'riST* l 


. -hrn 

fc,rj 2 

-hrn 

0 



Jl r i 2 T cl 

“ ^2 r il T c2 


(47) 


The r term comes directly from (10): 



Ti 3 

“ r i2 ' 


± 


2r,i 

*3 

II 

CO 

0 

0 

-7*il 


<P 

Q 

+ r„ 

~r t2 


. r *l 

0 




“7*13 


( 48 ) 
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and on combining all these, and substituting from the process noise model: 


H = 


(k 3 4 ^o)ri3 
0 

(lj 0 — k^)rn 


0 -2 U) 0 T iX 

(u; Q - ^l) r »3 -2u) 0 ri2 

(k\ 4 c^o) r il 0 


(fc 3 - 3r 0 )r i3 
0 

— (k 3 4 3r o )rii 


(3r o - k A )r i2 
(k4 4 3r o )rii 
0 


Cl 

C2 

C3 

rp 

9 



‘ -3r a 


+ r„ 

0 

+ k f 


r t3 

L 


*/^ 1 7*t 3 r C p2 + ^ 3 lr t2 r cpl 
-J 3 1 r<i^cpl “ J\^ r i3 r cp3 — m ~ 
J\ 1 »"i2 r cp3 ~ ^2 lr il r cp2 


[1 4 Wd(t)] 4 Vi 


(49) 


This completes the description of the accelerometers. There is 1 such 3 vector for each accelerometer. The 
noise depends critically on instrument design; but as we are interested only in feasibility, no particular 
instrument is used. Since a power spectrum is needed even for a generic instrument a cubic spectrum 
similar to (30) is assumed. The Appendix shows how the average power rt v (0), and the break frequency 
Wct) are determined from the rms acceleration error and the averaging time r of the measurement. 

5 Filter Structure 

The 1st step in calculating the terminal covariance in a dynamic estimation problem is to determine the 
structure of the filter. This starts with identifying the set of state variables that appear in the plant an 
measurement equations. From (22) and (25), it’s clear that we should choose: 

x = [ej, f 2 , «3> tfh <t>> 0\ T 

Following [7], it’s conventional to consolidate the plant equations in standard linearized form: 

x=Fx + G(u) + Bw ( 51 ) 

Here F is the plant matrix, u is a vector of controls, G(u), a possibly nonlinear vector function distributes 
the controls, w is a vector of independent process noises, and B is the process noise state distribution 
matrix. The matrices are readily identified. From (22) and (25), we find: 

ki 0 0 0 0 

0 0 0 fc 3 0 

0 0 0 0 k 4 

0 0 0 Wo 0 

1 0 -CJ 0 0 0 

010 0 0 

As for the control and process noise terms, it’s convenient to separate the deterministic process noise bias 
from the random components, and combine them with the actual controls, if any, to produce the G(u) 
used here. Since these terms will eventually cancel out in the analysis below the actual controls have no 
effect on filter performance, and there is no need to spell out G(u). Finally, by identifying w with w d { ) 
in (34), and including propeller torque, we have: 

B = kf[r cp $/Ji, r C p 2 / J 21 ~ r cpl/^3> 0, 0, 0] (53) 

Turning now to the measurement model, the direct appearance of the process noise in each of the ac- 
celerometer measurements requires a modification of the usual standard model: 

z = Hx + Yv + Uw + z B 

Here, H is the measurement partials matrix, developed above. From (49), this is: 

0 


F = 


0 

*2 

0 

1 

0 

0 


(52) 


(54) 


H,= 


(fc 3 + u} 0 )ri3 0 -2ui 0 rn 

0 (w 0 - fci)r i3 -2u> 0 r, 2 0 

(u> 0 — k 2 ) r il (£ 1 + u) o) r i2 0 


0 - 


(fc 3 - 3r o )r i3 (3r o - k 4 )r i2 

0 (k 4 + 3r o )ra 

(fc 3 + 0 


(55) 
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and the complete measurement partials matrix is: 

h=[ht, hi ...r (56) 

For example, if there are 7 vector accelerometers, H will be a 21 x 6 matrix. 

For measurement noise, it’s assumed that each axis of each accelerometer has separate independent noise 
Thus, v(t) has one element for each element of z, and Y is just an identity. A more elaborate model may 
be found in [1]; so Y is retained m what follows. The spectral properties of v(t) were developed above 
As for the process noise term, having established that w is w d (t), U comes immediately from (49): 


Ui = k f 


J 2 lr t3^cp2 + ^3 1 ^*t 2 ^*cpl 

“ ^3 r H r cpl - 1 Tt3T , cp3 — m _1 
^1 r i2 r cp3 “ J 2 * r il r cp2 


(57) 


cZr! V 1* K° 1Umn T t0r ? ith 3 SUCh elementS f ° r each Aerometer. The remaining terms in 
(49) constitute the bias z B . As it doesn’t affect the covariance analysis below, it’s not spelled out. 

An observer based on these models starts with an estimate x of the state x. This is caused to follow the 
deterministic parts of the plant equations (51), corrected by feeding back the residuals, i.e., the actual 
measurements z minus the measurement model (54). The filter structure then takes the form- 


* = Fx + G(u) + K(z - Hx - z B ) 


(58) 


Note that this structure assumes that the control and bias terms are known, and available to the filter 

dir " led here 13 l r hat G(u) 18 aCCUrately modeled i and that the biases have been accurately 
determined by some sort of in flight calibration. Pursuing these points is beyond our scope. 

6 Terminal Covariance 

rr:~^r amie filter is genera " > ' examincd by de,enninins ,he statis,ics ° f th ' *■ 

* = * ~ x (59) 

The evolution of £ comes from subtracting (51) from (58): 


(60) 


(61) 


£ = Z£ + KYv(t) - Ww(<) 

where the observer system matrix and the process noise effect matrix are defined by: 

Z = F-KH ; W = B - KU 

There’s lots to learn from (60). 1st, x, x, and all the control and bias terms have disappeared. Thus 

the “Wr r l eS f d0€ t n t dePend ° n . the controls - even if ^ey fail to stabilize the plant - 
all n ^ ^ eorem in the controls business. 2nd, filter stability requires Z to be stable- i e 

^ its eigenvalues are m the left half plane, a standard requirement in any negative feedback sy tern 

‘S differenUy: ^ 3 K Can bC f ° Und SUCh that Z is StabIe ’ ^en the state x fs saTd 
. e y the meas urements z. 3rd, the diagonal elements and the eigenvalues of Z have the 

eZvlT T ■T"" T:\ a " d filter ,ime i! “""“T ^ invL of i„ !ejt n^.i™ 

eigenvalue. This is used below to insure that the “optimal” filter has a reasonable settling time Finally 
since the noises are unbiased, so is £(t). 6 finally, 

Various measures have been proposed to study the quality of the estimate. Here, and generally in the 
references, attention has centered on the covariance of the error: 


e(0 = £[£(0£ t (01 


(62) 


where E is the expectation operator. The idea that, in a stationary system, P,(t) would have a terminal 
or asymptotic value, has been around a long time, but finding it could be quite tedious, if the settling time 
was ong. About 4 years ago, William McEneaney, in unpublished notes, showed that this terminaUalue 
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p € could be calculated directly from the structural information and the noise statistics. On generalizing 
to arbitrary power spectra, his ideas led to [6] and [7]. 

initial conditions have settled out. Then (60) is solved for “now in this form: 

(63) 


£(0) = r e Z "[KY v(m) - Bw((i)]d/i 
Jo 


t? <»r.ht :ix sl - 

x.^Thus, all terms of the form /i^x have the same dimensions, and if the exponential is merely viewe 
as shorthand for the formal expansion, there are no dimensional difficulties. 

The terminal covariance may now be found by substituting this into (62). 

P * = /°° f 0 eZ *{ KYE \ V ^ yT ^ YTKT + WE ( w M wT (* / )l WT } eZ VdfldV ^ 

This supposes that the expectation and integration operators may be commuted, and uses the assumption 
that w and v are independent and free of bias. On recognizing the autocorrelations of the noises, this . 

P C = r p e Z "[KYIU(M - u)Y t K t + WlUfc. - ( 65 > 

Jo Jo 

Well, autocorrelations and power spectra are Fourier transforms of each other. Using the one sided 
spectra of [6], these relations for any noise component are: 

rco 

fj( r ) = i [ S (u))c(Tw)du) ; S(u>)=2 I R(r])c(uT)dT (66) 

it Jo Jo 

After using the former in (65), and interchanging the order of integrations, there follows: 

p _ I / f f e Zfl Q (w)e ZT,/ c[a;(p - v)]dndvdu> 

4 71 Jo Jo Jo 


in which. q( w ) = KYS„(w)Y T K T + WS U) (w)W T 

Considerable progress can now be made by a change of coordinates: 

0 = /! + !' ; T)= H-V 

the double integration region is now the quadrant surrounding the +6 axis, so 

F J° e Zv/2 J°° e Z0/2Q( CJ )e ZTd/2 dOe- ZTr,/2 c(u>T))dT) 

+ f°° e Zi?/2 f e ’Z‘ 6 / 2 Q( u >)e ZTe/2 d0e~ Z v/2 c(u)i i)drt 

Jo Jv J 

Now, it’s not hard to establish that 

J e Ze/2 Q (uj)e ZTe/2 d9 = 2e Z<V2 D(w)e ZT<,/2 + constant 


(67) 

(68) 
(69) 


do: 


(70) 


(71) 
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where D(cv) satisfies the Lyapunov equation: 


ZD(cv) + D(<v)Z r = Q (tv) 


Putting this into (70), setting 77 
considerable simplification results: 


(72) 


-r) in the 1st integral, and evaluating at the required limits, a 


P< ~ •/ [ DM l + 


e^’ ? c(cv?j)d»7D(cv) I do; 


when another analytic integral has surfaced: 

J e^ v c(uT))df) = -(Z + <v 2 Z -1 )' 


(tv) 


(73) 


(74) 


leading finally to: 


1 f°° 

P « = ~y o [N(tv) + N T (tv))dtv ; N(tv) = (Z +tv 2 Z -1 ) -1 D(tv) 


(75) 

It may be noted that this analysis would break down in several places but for Z being stable Once attain 
especially in (74), the dimensions may look flaky. However, letting u, represent the dimensions of x ’ 
it is readdy shown from the differential equations that the expressions Z tJ t, Zy/ tv, and tvZ" 1 all have 

ItS— ? ,Ui , By e T Si0n ’ the ijth e,ement of < 74 > has the “o Thi J s work has 
established the forward procedure. For a given K, Z and W are computed from (61) and (61). A set of tv 

QMkfaenTt t0 C ° V !i tHe re ?' 0n Wh , Cre any ° f the n ° ise spectra are nonzero, with reasonable density. 
Q(cv) ,s then determined over this set from (68). Each Q(<v) yields a corresponding D(tv) by solution of 

he Lyapunov equation (72), and a corresponding N(tv) from (75). is then found by integrating (75). 

7 Optimal Feedback Gains 

Having found how to compute from K, we still need to find the K that yields optimal filter perfor 

SoblL W thCTe V am2i at ™ eans ' Whi,e P < C ° ntains the necessary information, in this 6th order 

problem there are 21 independent matrix elements; so some sort of scalar measure of P f is needed The 

software used here is based on a performance index 9 , constructed from the weighted trice of P, 


9 — •Pfaa/C'e 


(76) 


In this technique, known as “Bryson weighting”, each C t is the “level of concern” for the error A For 
example, if x, were a position, the level of concern might be C, = 1 m . C, = 10 m would show less 
S h T* optimization to put less weight on the variance of fc. Note that the Bryson 
technique has the virtue that , is the sum of dimensionless terms - it doesn’t add apples and oranges. 

kneads C toTz n wkb ^ perforrnance index ~ fiIter settling time. If the K that minimizes 

ime of t fiw u f ( h u ° Ugh n6gatlVe) ei « envalue ’ then we may see from (60) that the settling 

761 nenal P PS eXCessively S0 ' T « avoid such a problem, a term may be added to 

JltZ 0) ,f J f V™": To s “ h< ” “ d » * his ' •>» behavior o( .he filter evolution 

equation (60). If A Q symbolizes the eigenvalues of Z, and <r a = *(A a ), then the filter response to initial 

seTtUne Le'/Tc.' s',7 T be T'n' 1 ' “ ‘ of " exponentiaily decaying modes, with individual 
settling times -l/a a . Since all n modes decay simultaneously, the overall settling time is: 


t a = niax[— 9?(A q )] -1 = - JmaxK(A a )] * 


(77) 

Now suppose we introduce a concern level C t in seconds for the settling time t s . Then the overall 
performance index may be taken as: overall 


9 = {Pied Cl) + t s /C t 


(78) 
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we hive picked the concerns C, and C tl and have a K^ s uch that Z ^able, 

rr e,e^r: t sr* — « s 

a ;r^rrr P ^trr ££ «£ — - *- - 

*" uncontrolled gonn^ %"2f’ to know them to, say, 1%, our level of concern 
S£ ^ To-«^/»; a nd th2 is 25 - the toncen, level. However, an — * s.nngent rate Attar 
requirement, would shrink the rate concern levels. 

There i. one serious loose end - the Z ho"p«t"h 

on Kalman theory. Suppose each noise componen ( ) P van ; s hes for good. This level is 

the same average power R{ 0), and with cuto requency w infinity Replacing all the noise 

5 - 7rfi(0)/fl. The white noise “equivalent” to S(w) has level 5 out to intmity. nep * 

wtii|mnentl with these “equivalents” pauses d^n^nntctioVb^w^n^^nd^Pr^On 

25S5JJSE& - m - <«). - — *• - K explicil ’ we h ‘ ve: 

(79) 


where: 


KHP { + P«H t K t - FP« - P<F r = KMK t - KV t - VK T + BS„B 
M = YS v Y t + US„U T ; V = BS W U T 


(80) 




KM = P«H T + V 


(81) 


While this can’t be used directly to eliminate either K or P« from (79), we need only assume that some 
es ^ measurement componen. to Insure that M is non^gular. Thus: 


K = (P 4 H t + V)M 


-1 


(82) 


which, except for the V term, is a staple of Kalman theory. When this is substituted back into (79), an 
equally well known algebraic Riccati equation emerges: 

(83) 


A + XP$ + P<X T = P^H t M 1 HP« = P<LP< 


where: 


A = B(S„ - S 


u U 7 'M _1 US tt )B T 


x = F - VM ’H 


(84) 


All this reduces to Kalman theory when the me»umnen., drf. d^end on w(r); be U V^O. In 

-fl-r — * - *• 

Riccati equation has many solutions; but it’s known that at most 1 yields 
This is quite a large optimization. For example, if the 

K has 72 elements, all of which must be to alleviate this, 

aggravated by poor .conditioning ,m _^ or levels. On the hypothesis that the variance 

as 1 C? . oonsider paling the state variables and time: 


x'^Xi/Ci ; t' = t/C t 


(85) 
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lower case greek indices.' "T^vaVance of the s'ca^^b^L^then: 1118 P&Per ’ SUmmati ° n is ° nly over 

P £ij ~ E( x i x j) = P tijf{CiCj) (ggj 

*' = F'x' + G'(u) + B'w 

in which t° 7 l 

= FijCtCj/Ci ; B'ij = BijCt/Ci ; C' = G^/C, « 8 \ 

■* - * 


from which 


Hx = HV 

B '*i = CjHij 


The filter structure then becomes: 

* = F '*' + G'(u) + K'(z - H'x' - z B ) 
in which the derivatives are with respect to t', and 

Kj = 

and the error in the estimate 

(i = ~ x'i = £i/Ci 

evolves as 

i = zr + K'Yv(t) - W'w(t) 

where 

W' = B'-K'U ; Z' = F'-K'H' 

In components, these matrices are related to the unsealed versions by: 

Wl j =W ij C t /C i ; Zlj = ZijC t Cj/Ci 

AT ^ ^ 

foTe£::z m *:z by scaHng - «« the d — 


(89) 

(90) 

(91) 

(92) 

(93) 

(94) 

(95) 

(96) 


(97) 


=> a f a — C t a a t' = i*/C f 

On substituting these scaling relations into (78), , becomes rather simple; 

, = tV(P>J) + (i 

The modified iteration starts by forming B' and F' Then th* t r , 

y g and b ■ Then ’ the transformed algebraic Riccati equation is 

A' + X'P' +P 'X^=P'L'P' (99) 

in which A', X', and L' are computed as above except that F R «n/i w , , , 

equivalents. Note that V _ V' but no scaling «r X/f P * F ’ B ’ and H are replaced by their primed 

p; for the main iteration. Applying ^iSr * * — 


( 100 ) 


Q'H s K'YS„(o;)Y t K' t + W'S tu (ai)W ,T 

The Lyapunov equation is then: 

z'D'( w ) + D '( w )z' r = QV) non 

whose solution leads to N 7 fl nH P' u , 

unsealed values are V W " 9 haS Settled > yie,din S the terminal K' and P' , the 

P<ij=C i C j P{ ij ; Klj = Ci K-j/Ct 
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( 102 ) 


8 Results requires the exercise of several 

The calculation of the terminal covanance for a ^en s ^ ^ writt en in APL, and implemented 

programs in sequence, all more or less irrterac • P ^ h almogt entirely for the minimization 
on a 486DX 33 Mhz computer. A typical run requir AU the re8U u a cited here are based on 

of q, but including all the interactive input an ou p o n n 7 and 0 5 m; with a mass of 140 kg. A 
the numbers in the text. The spacecraft of 0.48458 rad/s. FYom 

momentum bias of 10 N-m^ » added, ^ g For numerical integration, 63 points were used in the 
the air data in the text, uw - 03604b r “/° r insure the acc uracy. The instrument 

« vector; but a couple of runs were repeated with more P«n£^ ^ drcumscribed sphe re has a 
consists of 4 accelerometers at the corners o gu 10 10 -s m / s 2 The averaging time was 1 s, for 
radius of 0.25 m. The noise levels ranged from 2 x 10 / 

= 62.832 rad/s. The C t are as in the text; C t - H» s. 

In all cases, q was dominated by t,; although kinematic equations, at 

Our interpretation is that * Ts to nofee increases, so must P„ and t. 

a natural frequency w OI and thus y a filter simulation would show that roll and pitch 

“ tab " 


/ 2 1 

10 9 x rms error - m/s 

0.2 

0.5 

1 

2 

5 

10 

a yaw - /xrad 
a roll - /irad 
g pitch - £trad 
is - S 

56.3 

3.04 

4.18 

223 

11.9 

3.33 

2.92 

315 

34.5 

7.06 

3.74 

336 

51.0 

32.0 
16.7 
670 

22.8 

20.1 

34.5 

695 

82.1 

72.5 

138 

1378 


_ „ rr otir We believe that this is due in part to the t s 
The progression to higher noiseseems ^ 3rd ^ comes from complex twins. The 1st 

dominance, but much more to th ( )• runs yielded quadruplets, composed 

and 4th run produced triplets, 1 of which iteration is unproductive. In all cases 

sss 

AoDendix — Averaged Measurement Noise , . . 

The instruments studied here are modeled “ • digit*' 

practice, they generally average the ana og ou p ments ^ tto! latcs T = 0. On the other hand, 

after each interval. The study ^ ( hdr d< , vtce ,',s delivering “samples” (really averages) 

the instrument manufacturers often characte R The noise assoc iated with these averages is 

every r seconds, or dt».fcly, « • ^^als with „,a,,„g this type of specification to 
then specified by a standard deviation a. Th ® a P p * situat ion was examined in [11], where it was 

the parameters of the assumed cubic power spectrum^ averages is; 

found that for an arbitrary noise power spectrum S(u), the variance 

<r 2 = 5 (w)[1-c(t«)]^5 (103) 

Jo 

Assuming the cubic spectrum (30) for the analog noise, the variance can be put in the form, 

*2 = R(0)/.(t« c ) 1 

where, in terms of the sine integral function: 

/.(*) = §Si(2*) + ^ (f + “) " ? (1 + S * X) 


(105) 


i 
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in which: 


Si (»)= f V — dz = y — 1 — 

Jo * 3 3! 5-5! 

The function looks ghastly for x < 1; but it actually behaves quite well: 


fs{x) = 1 - y +0(x 4 ) 


(106) 


(107) 


[MnvrH °Tr? I l g hm,t; ,e '’ ' f 8 time SeFieS is Veiy fre< l uentl y measured, but is long enough to cover 
many cycles of the highest noise frequency, then i?(0) is the variance of the samples, and the distinction 

isdsrcierli^^ aver /2 ge T?rr A / ctua ' iy ’ tws Hmit h ° ids f ° r ^ other i im it, , » i 

behavior can be sirTfrom the tabled * *' Vera ”’ A(x) ' S * monotonic decreasing function, whose 


X 

0 

0.1 

0.2 

0.5 

1 

2 

3 

fs(x) 

1 

0.99956 

0.99823 

0.98901 

0.9574 

0.84917 

0.71822 

X 

r / \ 

p 5 

^ 10 

n~ 20 

50 

100 

200 

r 500 

fs{x) 

0.50907 

0.28422 

0.14958 

.061632 

.031116 

.015633 

.006271 


When <7 was measured by the manufacturer, the repetition frequency 1/r was probably chosen about 
der of magnitude below the half power frequency c c /(2tt). Adapting this reasoning, we can pick: 


an 


so that rw c 20* = 62.832 rad; and R( 0) = .0492401a 2 . This assumed structure has been used to 
determine the measurement noise power spectrum in the study. 
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